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ABSTRACT 

Numerical experiments are performed to compare the accuracy obtained when 
physical and transform space filters are used to smooth the oscillations in Fourier 
collocation approximations to discontinuous solutions of a linear wave equation. 
High order accuracy can be obtained away from a discontinuity but the order is 
strongly filter dependent. Polynomial order accuracy is demonstrated when smooth 
high order Fourier filters are used. Spectral accuracy is obtained with the physical 
space filter of Gottlieb and Tadmor. 
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1. Introduction 


Spectral solutions to partial differential equations possess what is known as "spectral 
accuracy", that is, the error decays at a rate which depends only on the smoothness of the solution. 
If the solution is C°° on the domain then the error decays faster than any polynomial in the number 
of modes used (See Gottlieb, Hussaini and Orszag[2]). If the solution is analytic in a region about 
the real axis (which depends on the type of spectral expansion used), then the error decays 
exponentially [10]. In practical problems, error convergence which is faster than polynomial order 
is observed. Such behavior is seen, for example, when solving (smooth) two-dimensional 
nonlinear gas-dynamic problems[7]. 

If the exact solution of an equation has discontinuities, as is typical for hyperbolic PDEs, 
then the spectral solution is characterised by strong global oscillations. As a result, error 
convergence is poor even in smooth regions away from any discontinuities. The problem is 
compounded when solving nonlinear equations because these oscillations can lead to instabilities. 
The presence of global oscillations seems to indicate that spectral methods are inappropriate for the 
solution of hyperbolic equations when the solutions are discontinuous. 

Global oscillations in the discontinuous solutions can be removed by filtering. The main 
issue is whether or not "spectral accuracy" can be restored away from the discontinuity. In other 
words, can filtering be performed so that the error depends only on the local smoothness of the 
solution? 

The papers by Majda, McDonough and Osher[8] and Gottlieb and Tadmor[4] prove that for 
linear problems, proper smoothing allows spectral accuracy to be recovered in any region which 
does not include a discontinuity. The two approaches are different. Majda et al examined 
smoothing of Fourier approximations in frequency space. As a function of frequency, their filters 
were required to be unity in a neighborhood of zero frequency and smoothly decrease to zero at 
high frequencies. Gottlieb and Tadmor proposed smoothing in physical space with a highly 
oscillatory kernel whose support varies with position. Their procedure is applicable to spectral 
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methods based on other polynomial expansions such as Chebyshev and Legendre methods. Both 
papers include numerical experiments which show that with filtering, the solutions away from a 
discontinuity are far better than without. Neither shows, however, that the error is in fact spectral. 

There have been many attempts to apply filtering to the solution of nonlinear problems. 
Solutions of the Euler equations of gas-dynamics with shocks were presented in 
[1]>[6],[3],[9],[11],[13]. These papers indicate that some sort of filtering or artificial viscosity 
does give solutions which "look good". However, no one has presented solutions for which the 
error decay is spectral. Hussaini et al [6] tested several filtering strategies on a periodic problem 
with a shock wave and complex flow structure. They found that the best global error decay away 
from the shock was only first order. Solutions computed with the upwind finite difference scheme 
of Van Leer [12] were better than the spectral solution. 

In light of the fact that the papers by Majda et al and Gottlieb and Tadmor do not actually 
demonstrate spectral accuracy in practice and that the filtered solutions of nonlinear problems have 
systematically failed to exhibit spectral accuracy, we present numerical experiments on filtering. In 
particular we want to assess the usefulness of filtering for a realistic (rather than asymptotically 
large) number of grid points. We limit this paper to the study of Fourier collocation 
approximations to a linear model wave equation in order to test the effectiveness of filtering 
strategies for which theoretical evidence is available. We present results for both physical and 
transform space filters. 

2. Filtering of Fourier Spectral Approximations 

The Fourier collocation (also called pseudospectral) approximation to the linear model 
equation 
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with periodic boundary conditions u(0,t) = u(2tc,t) is constructed by approximating the derivatives 
u x at a number of grid (collocation) points. A uniform grid xj = 2icj/N is specified and 
approximate solution values vj are assigned to these points. A Fourier interpolant is passed 
through the N values of vj 
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The derivative approximations are computed by evaluating v'(xj) and we require that the 
differential equation is satisfied by vj at the grid points 
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Because these equations have constant coefficients, they can be solved exactly in time, so that any 
errors are due solely to the spatial approximation. 



(5) 


The initial condition for the differential initial value problem is u(x,o) = F(x). For the 
approximation we use vj = F(xj). This is appropriate for the scalar eq. 1 but not for systems [8], 

Fourier space filtering of the approximation for eq. 5 was proposed by Majda et al [8]. 
Here, the coefficients are modified by multiplying them by a filter function p(27tk/N). For the 
linear, constant coefficient problem, the filtered solution can be written 


( 6 ) 

v* ik(Xj-t) 

v j«) = -2/ ka kPk e 
-N/2 

The filter function p(tp) is a positive C°° cutoff function which vanishes near cp = ±ji. In a 
neighborhood of the origin, (p = 0, p = 1. With such filtering, Majda et a/[8] show that the error 
decay is spectral; that is, in any region R which does not include the discontinuity there is a 
constant C (which depends on X and the distance from a discontinuity) such that 

sup I u(x) - v(x) I < CN - ^ 

R 

for all X. 

Gottlieb and Tadmor[4] chose to filter in physical space. The rationale was that the Fourier 
space filtering, which is localized in transform space, includes the contribution of all points in 
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physical space, including points near a discontinuity. By filtering in physical space, the support 
of the filter can be defined so as to avoid these points. Smoothing is performed by convolving the 
collocation approximation with a regularization kernel 


(7) 

2n 

v*'¥ = J v(y) x P 0,p (x,y)dy 
o 

The kernel, \\f, is formed by the product of two factors. One is the Dirichlet kernel 

D p (y) = (1/2 tc) sin((p + l/2)y)/sin(y/2) (8) 

which is the physical space representation of truncation of the infinite Fourier series to a finite one 
of order p. The other factor localizes this in physical space. D p (y) is multiplied by a C 2s cutoff 

function, p, which satisfies p(0) = 1 and vanishes outside of the interval [-tc,tc]. Finally, a free 
parameter 0 is introduced so that the support of the kernel can be controlled. The resulting kernel 
is 


ve.P(y) = 0-lp(y/e)Dp(y/8) (9) 

For the collocation approximation, the convolution integral is replaced by the trapezoidal rule and 
the smoothed solution as a function of x is 
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Because the kernel is localized in physical space, terms for which (xj - y n )/0 is large will 
not contribute to the sum. Near the boundaries, i.e., xj near 0 and 2k, the sum must be modified 
to take into account the fact that the solution is actually periodic. Thus, (xj - y n -N)/0 may be 
within the support of \| f. To include points within the support of \|/ but not explicitly included in 
eq.10, the solution must be "wrapped around" to obtain accuracy near the boundaries. 

The error obtained by this smoothing is spectral. If p = (N/2)P where 0 < P < 1 then the 
error is bounded by 
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In other words, the error depends on the smoothness of the cutoff function. 

It is interesting to note that although the Dirichlet kernel has unit area on the interval, the 
regularization kernel does not. Thus, this smoothing does not preserve the mean of the function. 
If viewed in Fourier space, this means that the filter does not have amplitude identically one near 
the origin. As a consequence, even though the filter will asymptotically give spectral error decay, it 
can actually degrade the accuracy of a perfectly smooth function (such as a constant). A polynomial 
order Fourier space filtering will degrade the solution only in the higher modes, which will be quite 
small if the solution is smooth, . 

Both filtering strategies provide a lot of freedom. First, the Fourier space cutoff filter is not 
immediately defined. Such C°° functions exist, but are difficult to construct. Consequently, 
Majda et a/[8] suggested the use of an exponential cutoff filter 
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P(<P) = / r #1 , .2ml , , 

\ exp[-a(|<p|-q> c ) J cp c <|(p|<7t 

This filter satisfies the requirement that p = 1 in a neighborhood of 0 and decreases to zero at ± k , 

but it is no longer smooth. Majda et al preferred the fourth order filter, m =2. There is still 
freedom in choosing the cutoff frequency cp c and the coefficient, a. 

The freedom in specifying the physical space filter comes through the specification of the 
cutoff function, p and the parameters (3 and 9. Gottlieb and Tadmor[4] suggest the use of the C°° 
cutoff function 


P® = 


{ 


exp(a£ 2 /(£ 2 -l)) 

0 


1^1 < 1 
otherwise 


(13) 


The coefficient a is a free parameter. The parameter 0 determines the support of the filter. Gottlieb 
and Tadmor suggest a variable sized support which avoids the discontinuity. For a discontinuity at 
the point s, 0 = I x - si. 

3. Results of Numerical Experiments 

The effectiveness of physical and Fourier space filtering was tested on the model problem 


with the initial condition given by 
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0 < x < 2n 


which has a jump discontinuity at x = 0. An example of the solution when the discontinuity has 
moved to x = 7t and the sharpened cosine filter (see below) has been applied is shown in Fig. 1. 
The circles represent the spectral solution with 128 points and the solid line is the exact solution. 

It is important that the solution have high spectral content away from the discontinuity. If 
not, a filter may give results which are better than can be expected in general. For example, a first 
order filter can give excellent results away from a discontinuity if the solution is piecewise 
constant. 

Four Fourier space filters were applied to the solution of the linear model problem. Graphs 

of the amplitude vs. phase of these filters are shown on Fig. 2. The graph of the exponential 
cutoff filter, eq. 12, is for tp c = 7t/2. The coefficient was chosen so that the highest mode was 

multiplied by 10' 14, a number on the order of the machine epsilon. For comparison purposes, the 
use of the Lanczos filter in eq. 6 


p(cp) = sin(cp)/cp (15) 

is equivalent to differentiating in physical space with the second order centered difference 
approximation 

du/dx * (uj+i - uj_i)/2Ax (16) 

Another second order filter is the raised cosine filter 

p c (<p) = (1 + cos(<p))/2 (17) 
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This filter is equivalent to the common physical space average 


(uj+i + 2uj + uj_i)/4 

It can be sharpened to higher order by standard techniques [5]. An eighth order sharpened raised 
cosine filter is 


p = p c 4 (35 - 84p c + 70p c 2 - 20p c 3) (18) 

Results for Fourier space filtering are presented in figures 3-5. These are graphs of the 
logarithm of the absolute value of the pointwise error when the discontinuity was near x = k. For 

each filter, errors are presented for N = 32(circles), N=64(squares) and N = 128(diamonds). To 
keep the errors in perspective, 32 points are required to obtain errors on the order of 10‘ 5 for the 
closely related smooth problem with F(x) = (1 + cos 2 (6tcx/4))~1. 

Figures 3 and 4 are graphs of errors smoothed with the Lanczos and sharpened raised 
cosine filter. As expected, the error near the discontinuity is 0(1), regardless of the filter used. 
Since these are semi-log plots, a uniform spacing between the curves for different values of N 
indicates polynomial accuracy. Both filters give the expected results: The spacing between the 
curves for the Lanczos filtering is roughly 0.6 which corresponds to second order convergence. 
Eighth order accuracy is observed for the sharpened raised cosine filter. 

Under the same circumstances, the exponential cutoff does not give spectral accuracy. Fig. 
5 shows the results of the fourth order filter with a cutoff at 0.5 tl Though the error appears to 
decrease faster than polynomial order in some areas, the error is not smooth. Furthermore, the 

region of large error about the jump is also quite large. If all the modes except k = 0 are smoothed 
((p c = 0) then the region of high error about the jump shrinks to that of the other polynomial order 

filters. The error decay, however, is strictly fourth order. These results are consistent with what 
was observed by Majda et al[ 8]. 


9 



Spectral accuracy has been obtained with physical space filtering. Parameters close to 
those proposed by Gottlieb and Tadmor[4] were used. The errors do not appear to be sensistive to 
the parameter a and a = 10 is illustrated. They are sensistive to |3 and the best results were 
obtained for (3 = .9. 

Figure 6 is a graph of the errors for the physical space filter. Away from x=7t the error 
decay ranges between fourth and eighth order as the number of modes increases and hence is 
spectral. Still, this filter is outperformed by the eighth order sharpened cosine filter. This is pardy 
due to the fourth order behavior for small N. But for fixed N the error varies as much as four 
orders of magnitude within the smooth region. Errors near x=k are 0(10) rather than 0(1). This 
behavior near the jump is consistent with the error estimate, eq. 12, which gives a large error 
bound when 0 is very small. In the smooth parts of the solution, the error is largest where the 
parameter 0 is largest. 

The quality of the physical space filter is very dependent on 0. For example, the solution 
can be substantially improved on this problem by n£l using variable support in x. Fig. 7 shows 
the errors for the constant value of 0 = 1 , which gives the best results for this problem. The errors 
near the jump are now 0(1). For each N the error is now more uniform over the smooth region 
and the error decay ranges between fourth and tenth order. Unfortunately, the improvement seen 
in the smooth part is problem dependent. For the initial condition F(x) = (1 + (x/rc)2)-l, the errors 
away from the jump are an order of magnitude better using the variable 0 than using a constant 0. 

Keeping the physical space support of the filter fixed in space does have an advantage in 
the efficiency of the filtering method. The convolution represented by eq. 10 requires O(N^) 
operations to compute. Represented in transform space, this can be done with a fast Fourier 
transform in O(NlogN) operations, but only if 0 is kept constant. 

The results presented here show that the asymptotic decay of the error is very filter 
dependent. For coarse grids, e.g. N = 32, however, the differences are not great. Figs. 3-7 
indicate that the error after any of the high order smoothings is only roughly a third of a second 
order finite difference method as modelled by the Lanczos filtering. The lack of clear superiority of 
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the spectral method over a low order finite difference method for coarse meshes in the presence of 
discontinuities appears to be typical. 

4. Conclusions 

The results of this investigation show that it is possible to obtain high order solutions to the 
discontinuous linear model initial value problem. Even though the Fourier approximation is first 
order everywhere, away from a jump discontinuity it is possible to extract higher order 
information by filtering. The order of accuracy is very filter dependent. Smooth, high order 
polynomial filters, of which eq. 18 is an example, are easy to construct and their behavior is 
typical. Sacrificing smoothness for infinite order tangency at zero frequency as proposed by Majda 
et al does not provide satisfactory results. Spectral accuracy has been obtained with the filter of 
Gottlieb and Tadmor in the sense that the order of accuracy of the method increases as the number 
of points increases. For a coarse mesh, however, there is little difference between a second order 
finite difference method and a computationally more expensive spectral method. 

Of the filters tested, the spectrally accurate physical space filter of Gottlieb and Tadmor 
gives the worst results near a discontinuity. This could be disastrous for a nonlinear problem. The 
performance of the filter can be improved by allowing the the support of the filter to include the 
discontinuity. For the problem described here, the best overall results away form the jump were 
obtained by keeping 9 = 1 everywhere, but this is not always the case. 

Finally, we remark that the conclusions of this paper do not necessarily carry directly over 
to the solution of nonlinear problems. Though eighth order accuracy is observed for the filter eq. 
18, Hussaini et al[6] observed only first order accuracy when it was applied to the solution of a 
nonlinear problem. 
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FIGURE CAPTIONS 


1. Graph of the solution of eq. 1 with initial condition, eq. 14. The solid line is the exact solution, 
circles are computed solutions smoothed with the sharpened raised cosine filter. 

2. Fourier space filters, a: raised cosine; b: Lanczos; c: sharpened raised cosine; d: exponential 
cutoff. 

3. Pointwise errors for the solution filtered with the second order Lanczos filter. 

4. Pointwise errors for the solution filtered with the eighth order sharpened raised cosine filter. 

5. Pointwise errors for the solution filtered with the exponential cutoff filter. 

6. Pointwise errors for the solution filtered with the variable support physical space filter. 

7. Pointwise errors for the solution filtered with a constant support physical space filter. 
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